#find x to maximized f
#f = (1+bx)^p(1-x)^q
#ifelse(g,bx,-x),P(g==T) = p, P(g==F) = q

n = 100
p = 0.7 # 
b = 0.8 #return from investing 1 unit, odds b to 1
kf = (b*p-(1-p))/b #Kelly formula, the fraction of total capital
if(kf < 0){
   kf = 0
   warning("warning: not to invest ")
}
x = 0.02
er = 1+x*(p*b+p-1)
#(1+b*x)*p+(1-x)*(1-p) 
#p+b*p*x+1-x-p+x*p = 1-x+p*x*b+p*x

vr2 = x^2*(1+b)^2*p*(1-p)
#(1+b*x - 1+x-p*x*b-p*x)^2 * p +(1-x - 1+x-p*x*b-p*x)^2*(1-p) 
#((1-p)*x*(b+1))^2*p + (p*x*(1+b))^2*(1-p)

sdr = x*(1+b)*sqrt(p*(1-p)) # linear to x

assumed you invest 20% of your capital. then after one period, your expected capital is 1.052 with standard deviation 0.1649727

#the function kelly try to optimized -- the accumulated return.
# as the number of bets increase, W/N -> p, and L/N -> 1-p. so it converge to this function
zf = function(r,b=1.6,p=0.6){(1+b*r)^p*(1-r)^(1-p)}


#in term of the tree method, for n bets
zf.tree = function(x,b=1.6,p=0.6,n=100){
   ss = c(1+b*x,1-x)
   ps = dbinom(0:n,n,p)
   A = c(1,rep(ss[1],n)) %>% cumprod()
   B = c(1,rep(ss[2],n)) %>% cumprod() %>% rev()
   return(list(ps = ps, values = A*B))
}


#in term of the tree method, for n bets
zft = function(x,b=1.6,p=0.6,n=100){
   ans = zf.tree(x,b,p,n)
   ps = ans$ps
   Vs = ans$values
   e.v = sum(ps*Vs)
   e.s2 = sum(ps*(Vs-e.v)^2)
   return(list(e.ret = e.v, e.sd = sqrt(e.s2)))
}

#in term of the tree method, for n bets
zft.log = function(x,b=1.6,p=0.6,n=100){
   ans = zf.tree(x,b,p,n)
   ps = ans$ps
   Vs = log(ans$values) #use log transform, because the upper end is too big. take it as utility function
   e.v = sum(ps*Vs)
   e.s2 = sum(ps*(Vs-e.v)^2)
   return(list(e.ret = e.v, e.sd = sqrt(e.s2)))
}

zft.mr = function(x,b=1.6,p=0.6,n=100,a=0.95){
   ans = zf.tree(x,b,p,n)
   psc = cumsum(ans$ps)
   e.v = ans$values[which.max(ans$ps)]
   a1 = (1-a)/2
   a2 = 1-a1
   ind1 = which.min(abs(psc -a1))
   ind2 = which.min(abs(psc -a2))
   e.sd = ans$value[ind2]-ans$value[ind1]
   return(list(e.ret = e.v, e.sd = e.sd, e.sd_low = ans$value[ind1], e.sd_hi = ans$value[ind2]))
}

zft.mm = function(x,b=1.6,p=0.6,n=100,w=0.5){
   ans = zf.tree(x,b,p,n)
   ind = which.max(ans$ps)
   e.v = ans$values[ind]
   c1 = 0
   ind1 = ind-1
   while((ind1>0) & (c1 < w/2)){
      c1 = c1+ans$ps[ind1]
      ind1 = ind1-1
   }
   c2 = 0
   ind2 = ind+1
   while((ind2<n+1) & (c2 < w/2)){
      c2 = c2+ans$ps[ind2]
      ind2 = ind2+1
   }
   if(ind1 == 0) ind1 = ind1+1
   if(ind2 >n) ind2 = ind2-1
   e.sd = ans$value[ind2]-ans$value[ind1]

   #print(c(ind,ans$ps[ind],ind1,ind2,ans$ps[ind1:ind2]))
   return(list(e.ret = e.v, e.sd = e.sd, e.sd_low = ans$value[ind1], e.sd_hi = ans$value[ind2]))
}

#x is the fraction to be used for the bet
gdx_data = function(p,b,x,m_paths=500,n_steps = 100){
   ps = c(p,1-p)
   ss = c(1+b*x,1-x)
   r.ind = sample(c(1,2),m_paths*n_steps,replace = T, prob = ps)
   r = ss[r.ind]
   gd = data.frame(path = gl(m_paths,n_steps),d.ret = r,o = rep(1:n_steps,m_paths))
   gd = gd %>% group_by(path) %>% mutate(c.ret = cumprod(d.ret))
   gdx = rbind(data.frame(path = gl(m_paths,1),
                          d.ret = rep(1,m_paths),
                          o = rep(0,m_paths),
                          c.ret = rep(1,m_paths)),gd)
   return(gdx)
}
data = gdx_data(p,b,kf,m_paths = 800)
data$x = kf
xs = c(0.02,0.2,0.5,0.7,0.8)
for(x in xs){
   gdx = gdx_data(p,b,x)
   gdx$x = x
   data = rbind(data,gdx)
}
gp = ggplot(data = data,aes(x = o, y = c.ret, group = path))  + 
      geom_line(alpha = 0.02) +scale_y_log10()+facet_wrap(~x)
print(gp)

the theoratical distribution of different fractions to invest

x = 0.2
gf = function(x){zf.tree(x,b,p) %>% as.data.frame() %>% mutate(x = x, pc = cumsum(ps),cf95 = pc > 0.025 & pc <0.975 )}
g = gf(x)
ps = g$ps
ind = which.min(abs(cumsum(ps)-0.975))
p.level = ps[ind]
xs = c(0.02,0.2,0.5,0.7,0.8,kf)
g.data = map_dfr(xs,gf)
pg = ggplot(g.data %>% filter(values > 1e-15 & values < 1e6),aes(x = values, y = ps,color = factor(x)))+
   geom_point(alpha = 0.7,size = 0.5)+geom_line(aes(group = factor(x)),alpha = 0.5)+
   geom_vline(xintercept = 1,linetype="dashed")+
   geom_hline(yintercept = p.level,color="red",alpha = 0.4,linetype="dashed")+
   #facet_wrap(~factor(x),ncol = 1)+
   scale_x_log10()+theme(legend.position = "none")
print(pg)

a few points

data.end = data %>% filter(o == max(o))  %>% mutate(log.c.ret = log(c.ret))
data.end.stat = data.end %>% group_by(x) %>% 
      summarise(e.ret = mean(c.ret),s.ret = sd(c.ret),
                m.ret = quantile(c.ret,0.5),
                q05.ret = quantile(c.ret,0.05),
                q95.ret = quantile(c.ret,0.95),
                e.elog.ret = mean(log.c.ret)%>%exp(),s.elog.ret = sd(log.c.ret)%>%exp(),
                .groups = "keep")
knitr::kable(data.end.stat,digits = 3)
x e.ret s.ret m.ret q05.ret q95.ret e.elog.ret s.elog.ret
0.020 1.679 0.278 1.657 1.242 2.137 1.656 1.180
0.200 129.715 288.974 40.249 2.940 786.508 36.624 5.153
0.325 1491.265 6923.664 80.362 0.545 6346.273 77.107 15.409
0.500 21932.848 122078.999 15.778 0.011 59611.065 21.489 121.295
0.700 2614.954 32623.556 0.007 0.000 698.848 0.003 2268.210
0.800 666.551 8393.498 0.000 0.000 2.929 0.000 20291.095
gp = data.end %>%
      ggplot(aes(x = x, y = c.ret))  + 
      stat_summary(alpha = 0.8,geom = "bar",fun = "mean",width = 0.03)+
      stat_summary(geom = "errorbar",fun.data = "mean_se",width = 0.01)+
      geom_hline(yintercept = 1,linetype = "dashed")+
      scale_y_log10()
print(gp)

#p = 0.7 # 
#b = 0.7 #return from investing 1 unit, odds b to 1
#er = b*p-1*(1-p) #expected return
#kf = (b*p-(1-p))/b 

xs = c(0.02,0.2,0.5,0.7,0.8)
zx = seq(0,0.95,0.01)
zd = data.frame(fraction = zx, expected.return = zf(zx,b,p),sdr = zx*(1+b)*sqrt(p*(1-p))/10)
zd2 = data.frame(fraction = xs, expected.return = zf(xs,b,p),sdr = xs*(1+b)*sqrt(p*(1-p))/10)
zd$color = "grey"
zd2$color = "blue"
zd = rbind(zd,zd2)

pg = ggplot(data = zd)+geom_point(aes(x = fraction, y = expected.return,color = color))+
         geom_line(aes(x = fraction, y = expected.return+sdr),color = "grey")+
         geom_line(aes(x = fraction, y = expected.return-sdr),color = "grey")+
         geom_point(aes(x = kf,y = zf(kf,b,p)),color = "red")+
         scale_colour_identity()
print(pg)

#p = 0.7 # 
#b = 0.7 #return from investing 1 unit, odds b to 1
#er = b*p-1*(1-p) #expected return
#kf = (b*p-(1-p))/b 

xs = c(0.02,0.2,0.5,0.7,0.8)
zx = seq(0,0.8,0.01)
f.zft = function(x){zft.log(x,b,p)  }

zd =  map_df(zx,f.zft) %>% rename(expected.return = e.ret, sdr = e.sd) %>% 
         mutate(fraction = zx,sdr = sdr,color = "grey")
zd2 =  map_df(xs,f.zft) %>% rename(expected.return = e.ret, sdr = e.sd) %>% 
         mutate(fraction = xs,sdr = sdr,color = "blue")

zd = rbind(zd,zd2)

pg = ggplot(data = zd)+geom_point(aes(x = fraction, y = expected.return,color = color))+
         geom_line(aes(x = fraction, y = expected.return+sdr),color = "grey")+
         geom_line(aes(x = fraction, y = expected.return-sdr),color = "grey")+
         geom_point(aes(x = kf,y = f.zft(kf)$e.ret),color = "red")+
         scale_colour_identity()+ylab("expected.log.return")
print(pg)

pg = ggplot(data = zd)+geom_point(aes(x = fraction, y = expected.return %>% exp(),color = color))+
         geom_line(aes(x = fraction, y = (expected.return+sdr/10)%>% exp()),color = "grey")+
         geom_line(aes(x = fraction, y = (expected.return-sdr/10)%>% exp()),color = "grey")+
         geom_point(aes(x = kf,y = f.zft(kf)$e.ret %>% exp()),color = "red")+
         scale_colour_identity()+ylab("expected.log.return")
print(pg)

#p = 0.7 # 
#b = 0.7 #return from investing 1 unit, odds b to 1
#kf = (b*p-(1-p))/b 
#if(kf < 0) kf = 0
xs = c(0.02,0.2,0.5,0.7,0.8)
zx = seq(0,0.8,0.01)
f.zft = function(x){zft.mr(x,b,p,a = 0.3)} #the central 30% confident zone

zd =  map_df(zx,f.zft) %>% rename(expected.return = e.ret, sdr = e.sd) %>% 
         mutate(fraction = zx,sdr = log1p(sdr),color = "grey")
zd2 =  map_df(xs,f.zft) %>% rename(expected.return = e.ret, sdr = e.sd) %>% 
         mutate(fraction = xs,sdr = log1p(sdr),color = "blue")

zd = rbind(zd,zd2)

pg = ggplot(data = zd)+geom_point(aes(x = fraction, y = expected.return,color = color))+
         geom_line(aes(x = fraction, y = e.sd_low),color = "grey")+
         geom_line(aes(x = fraction, y = e.sd_hi),color = "grey")+
         geom_point(aes(x = kf,y = f.zft(kf)$e.ret),color = "red")+
         scale_colour_identity()+ylab("most.likely.return")


print(pg+scale_y_log10())

print(pg)

the relationship between variance and return

zx = seq(0.01,0.8,0.01)
f.zft = function(x){zft(x,b,p)}
d = map_df(zx,f.zft) %>% mutate(fraction = zx)
pg = ggplot(d)+geom_point(aes(x = e.sd,y = e.ret,color = fraction))+
      scale_y_log10()+scale_x_log10()
print(pg)

zx = seq(0.01,0.8,0.01)
f.zft = function(x){zft.log(x,b,p)}
d = map_df(zx,f.zft) %>% mutate(fraction = zx)
pg = ggplot(d)+geom_point(aes(x = e.sd,y = e.ret,color = fraction))+
      geom_point(aes(x = f.zft(kf)$e.sd, y = f.zft(kf)$e.ret),color = "red")

print(pg)

zx = seq(0.01,0.8,0.01)
f.zft = function(x){zft.mr(x,b,p,a = 0.95)}
d = map_df(zx,f.zft) %>% mutate(fraction = zx)


pg = ggplot(d)+geom_point(aes(x = e.sd,y = e.ret,color = fraction))+
      geom_point(aes(x = f.zft(kf)$e.sd, y = f.zft(kf)$e.ret),color = "red")

print(pg)

zx = seq(0.01,0.8,0.01)
f.zft = function(x){zft.mr(x,b,p,a = 0.3)}
d1 = map_df(zx,f.zft) %>% mutate(fraction = zx)
d1$a = 0.3

f.zft = function(x){zft.mr(x,b,p,a = 0.6)}
d2 = map_df(zx,f.zft) %>% mutate(fraction = zx)
d2$a = 0.6

f.zft = function(x){zft.mr(x,b,p,a = 0.1)}
d3 = map_df(zx,f.zft) %>% mutate(fraction = zx)
d3$a = 0.1

d = rbind(d1,d2,d3)

pg = ggplot(d)+geom_point(aes(x = e.sd,y = e.ret,color = fraction))+
      geom_path(aes(x = e.sd,y = e.ret,linetype = factor(a)),color = "grey")

print(pg)

zx = seq(0.01,0.8,0.01)
f.zft = function(x){zft.mm(x,b,p,w = 0.9)}
d = map_df(zx,f.zft) %>% mutate(fraction = zx)


pg = ggplot(d)+geom_point(aes(x = e.sd,y = e.ret,color = fraction))+
      geom_point(aes(x = f.zft(kf)$e.sd, y = f.zft(kf)$e.ret),color = "red")

print(pg)

zx = seq(0.01,0.8,0.01)
f.zft = function(x){zft.mm(x,b,p,w = 0.3)}
d1 = map_df(zx,f.zft) %>% mutate(fraction = zx)
d1$w = 0.3

f.zft = function(x){zft.mm(x,b,p,w = 0.6)}
d2 = map_df(zx,f.zft) %>% mutate(fraction = zx)
d2$w = 0.6

f.zft = function(x){zft.mm(x,b,p,w = 0.1)}
d3 = map_df(zx,f.zft) %>% mutate(fraction = zx)
d3$w = 0.1

d = rbind(d1,d2,d3)

pg = ggplot(d)+geom_point(aes(x = e.sd,y = e.ret,color = fraction))+
      geom_path(aes(x = e.sd,y = e.ret,linetype = factor(w)),color = "grey")

print(pg)

the fair bet

DiagrammeR::grViz("digraph {
  node [shape = oval]
  a [label = '1']
  b [label = '1+bx']
  c [label = '1-x']
  a->b [label='p']
  a->c [label='1-p']
  
  d [label = 'x']
  e [label = '(1+b)x']
  f [label = '0']
  d->e [label='p']
  d->f [label='1-p']
  {rank=same;a;d;}
  }",
  height = 200)
N = 1e3
#pb.data = data.frame(p = runif(N),b = rnorm(N,mean = 1,sd = 0.3)) %>% filter(b >0)
pb.data = data.frame(p = runif(N),b = runif(N,min = 0,max = 2))%>% filter(b >0)

pb.data = pb.data %>% mutate(kf = (b*p-(1-p))/b,jd = kf > 0,kf_ = ifelse(kf>0,kf,0))
pg = ggplot(pb.data)+geom_point(aes(x = p, y = b,shape = jd,color = kf_))+
      geom_line(data = data.frame(x = seq(0.33,1,0.01)),aes(x = x,y = 1/x-1),color = "red")+
      guides(shape=guide_legend(title="kf > 0"))+
      scale_shape_manual(values=c(1,3))
print(pg)